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Abstract 


Although liquid water has been studied for many decades by (X-ray and neutron) diffraction measurements, new 
experimental results keep appearing, virtually every year. The reason for this is that neither X-ray, nor neutron diffraction 
data are trivial to correct and interpret for this essential substance. Since X-rays are somewhat insensitive to hydrogen, 
neutron diffraction with (most frequently, H/D) isotopic substitution is vital for investigating the most important feature 
in water: hydrogen bonding. Here, the two very recent sets of neutron diffraction data are considered, both exploiting the 
contrast between light and heavy hydrogen, and ^H, in different ways. Reverse Monte Carlo structural modeling is 
applied for constructing large structural models that are as consistent as possible with all experimental information, both 
in real and reciprocal space. The method has also proven to be useful for revealing where possible small inconsistencies 
appear during primary data processing: for one neutron data set, it is the molecular geometry that may not be maintained 
within reasonable limits, whereas for the other set, it is one of the (composite) radial distribution functions that cannot 
be modeled at the same (high) level as the other three functions. Nevertheless, details of the local structure around the 
hydrogen bonds appear very much the same for both data sets: the most probable hydrogen bond angle is straight, and 
the nearest oxygen neighbours of a central oxygen atom occupy approximately tetrahedral positions. 
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1. Introduction 

Liquid water, as it is the the basis of life on Earth, 
is a most common, everyday liquid. It is the simplest 
compound of two universal elements, hydrogen and oxy¬ 
gen (H 2 O), and is the second most frequent molecule in 
the Universe. In spite of these, water is still an amazing 
substance with unique properties: it has phas^ density, 
thermodynamic and other physical anomalies [l|. Hence 
it is not surprising that water has been, and still is the 
most researched liquid. (For a comprehensive summary, 
see e.g. the collection of references in [l|.) 

In order to understand its properties (and anomalies!), 
it is necessary to comprehend the microscopic structure of 
liquid water. There are numerous studies in the literature 
using different spectroscopy and scattering techniques: X- 
ray [3, ill and neutron scatte ring j i i i , X-ray ab¬ 


sorption spectroscopy (XAS) 
spectroscopy (XES) [iSl UA 


I , X-ray emission 
IH I or small-angle X-ray 


scattering experiments (SAXS) [17|,118j. To interpret the 
experimental results many models and computer simula¬ 
tion works have been performed (see e.g. iHiSlIllil 
and references therein). 
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Though there have been many investigations, there is 
not a reassuring consensus present about the structure of 
water. For example, even the concept of tetrahedral ar¬ 
rangement of 4-fold coordinated water molecules, that had 
been widely accepted earlier, has been challenged and an 
asymmetric distribution of twofold-coordinated molecules 
was suggested [3|. The proposed structural organization, 
with hydrogen-bonded chains and rings of water molecules 
in a weakly hydrogen-bonded disordered network, is still 
controversial (see e.g. 13, 24|). As an interesting addition 
to the ongoing debate concerning coordination numbers, 
Skinner et al. report an ’isosbestic point’, a distance at 
which the integrated 0-0 coordination number is inde¬ 
pendent of temperature!^. 

Diffraction measurements are the obvious means for 
structure determination. X-ray diffraction can be useful 
for determining the oxygen-ox yge n (and, to some extent, 
oxygen-hydrogen) correlations [23, but this technique is 
not too sensitive for hydrogen. Since neutron scattering 
lengths are different for the hydrogen isotopes, neutron 
diffraction with H/D isotope substitution is, in principle, 
suitable for the separation of the hydrogen-hydrogen and 
hydrogen-oxygen correlations (if one accepts the approxi¬ 
mation that the structure of water depends only weekly on 
the isotopic composition). In practice, however, the cor¬ 
rection of the measured scattering data is extremely diffi- 
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cult due to the strong incoherent and inelastic scattering 
of neutrons by protons (^H). There is still not a generally 
applicable way to remove such effects from the measured 
data (see, e.g., This is why structural model¬ 
ing [ 2 ^ has played an important role in the 

history of structure studies on liquid water. 

One of the possibilities that can be useful for interpret¬ 
ing the corrected diffraction data sets is to prepare three- 
dimensional atomic (structural) models that are consistent 
with all of the input data sets. The Reverse Monte Carlo 
(RMC) method is a perfect tool for this purpose. RMC 
also allows one to check whether it is possible to prepare 
physically meaningful particle confi gura tions that match 
corrected experimental data (see Ref.[29|); i.e., the method 
(under certain conditions) may serve for providing a kind 
of ’first aid support’ to measured diffraction data. 

In this work we have studied two sets of very recent 
neutron total scattering structure factors (TSSF) and the 
corresponding radial distribution functions. One of them 
0, referred to as ’Case 1’ throughout, presented diffrac¬ 
tion measurements on four mixtures of light and heavy wa¬ 
ter; the partial radial distribution functions (PRDF) were 
also provided. The other study Q, indicated as ’Case 
2’ in the following, is based on neutron diffraction mea¬ 
surements with oxygen isotope substitution. From these 
measurements the authors obtained linear combinations of 
the 0-0 and 0-H (or 0-D) partial structure factors. They 
presented radial distribution functions obtained by Fourier 
transformations, as well. 

Two important differences between the two sets of neu¬ 
tron diffraction experiments should be mentioned explicitely: 

(1) Data for ’Case 1’ have been gathered at a spallation 
neutron source (ISIS, UK), whereas experiments connected 
to ’Case 2’ were performed on a steady-state reactor based 
instrument (ILL, France); this difference has a profound 
effect on the methods and complexity of data treatment 
(see Refs. S0)- 

(2) Since it was not stated otherwise in Ref.Q, stan¬ 
dard values of the atomic coherent scattering lengths were 
assumed for data in ’Case 1’, whereas for ’Case 2’, they 
were taken from the original publication^. The coherent 
scattering lengths applied throughout this work are quoted 
in Table [TJ 



bnatQ 

bisQ 


bn 

Case 1 

5.803 

- 

6.671 

-3.7406 

Case 2 

5.805 

6.005 

6.619 

-3.7409 


Table 1: Coherent atomic scattering lengths used for neutron data 
of ’Case I’Q and ’Case 2’0| (in fm). 


We performed Reverse Monte Carlo modeling, in or¬ 
der to (a) learn about the internal consistency of these 
new, much improved scattering data and (b) find struc¬ 
tural models that would be consistent with the entire data 


sets, including the r-space information, simultaneously. 


2. Reverse Monte Carlo modeling 

Reverse Monte Carlo (RMC) modeling has been de¬ 


scribed in many publications in detail (30l l3ll . l32l l33| . 


so here only a brief description is necessary. RMC is an 
inverse method to obtain large three-dimensional struc¬ 
tural models that are consistent with the supplied (ex¬ 
perimental and/or theoretical) data sets (TSSF-s and/or 
(P)RDF-s). It can be used in conjunction with any quan¬ 
tity that can be expressed directly from the atomic coor¬ 
dinates. By moving particles randomly in the simulation 
box, the difference between the experimental (or ’quasi- 
experimental’, see, e.g., and RMC model structural 
quantities (e.g. structure factors) is minimized. As a re¬ 
sult, particle configurations are obtained that are consis¬ 
tent with all the input data. From the particle config¬ 
urations, further structural characteristics (coordination 
numbers, nearest neighbor distances, bond angle distribu¬ 
tions, etc...) may be calculated. Over the past nearly 
three decades, RMC has been successfully applied to a 
wide variety of systems, from simple liquids (see, e.g., the 
review of McGreevy (3flllL t hrough metallic [35| and cova¬ 
lent [ 3 ^ glasses, simple [s^, [s^ and hydrogen-bonded Hi 
molecular liquids, as well as for characterizing disordered 
crystalline structures (11,0. 

The method is also suitable for establishing whether 
various input data sets are consistent with each other: if 
they are then they can be fitted simultaneously within 
their uncertainties [3. If, on the other hand, not each 
element of the input set is consistent with the others then 
it is possible to tell which element is problematic: an ap¬ 
proach proposed and tested recently for input data 
consisting of TSSF-s and PRDF-s is applied in this work. 

Here the RMC-|--|- codejs^ was used to obtain struc¬ 
tural models. The cubic simulation boxes contained 6000 
atoms (2000 molecules), the atomic number density was 
0.1 A“^, the simulation box length was 39.1 A. In order to 
keep the atoms together in the molecules during the calcu¬ 
lations the ’fixed neighbor constraints’ (FNC) option was 
applied [^. This algorithm connects two hydrogen atoms 
and their oxygen central atom permanently via their iden¬ 
tity numbers (defined for the atomic configuration file). 
For the realization of this constraint, the 0-H and H-H 
intramolecular distances were kept between minimum and 
maximum values: the values for the various FNC combi¬ 
nations tested are shown in Tabled] Intermolecular closest 
approach (’cutoff’) distances between atoms were chosen 
as indicated in Table dl 

Simulations were started from a random configuration 
where only the cutoff distances and the fixed neighbor con¬ 
straints were in effect (’hard sphere’ model). The number 
of accepted moves was about 1 — 2 x 10^ in each calcu¬ 
lation. a parameters (essentially, control parameters for 
the different data sets that influence the tightness-of-fit. 
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0-0 

O-H 

H-H 

intramolecular 
distances 
realistic molecule 

- 

0.95 - 
1.02 

1.50 - 
1.62 

intramolecular 
distances ’elastic’ 


0.78 - 
1.15 

1.36 - 
1.71 

molecule 


intermolecular 




distances 
all models 

2.3 

1.5 

1.6 


Table 2: Minimum and maximum O-H and H-H intramolecular dis¬ 
tances used in fixed neighbor constraints and intermolecular closest 
approach distances (cutoff distances) between atoms (in A). 

cf. Ref. 0) were decreased progressively during the sim¬ 
ulation runs, resulting in a gradually improving fit to the 
target functions. Final a values are listed in Tables [3] and 
m Note that these tables also provide information on the 
sets of input data applied in the various RMC calcula¬ 
tions: for Case 1 (data from Ref. @), always at least 
7 target functions (4 TSSF-s, 3 intermolecular PRDF-s) 
were considered, whereas for Case 2 (data from Ref. 0), 
4 functions (two in Q and two in r space, each a kind of 
’composite’) were approached. 



cr 

R - 
Case 
la 

R - 
Case 
lb 

R - 
Case 

Ic 

gooir) 

0.02 

6.2 

6.1 

6.1 

goH{r) 

0.02 

5.9 

4.6 

5.1 

gHHir) 

0.02 

4.1 

3.0 

4.1 


0.02 

- 

- 

15.4 

Fo^oiQ) 

0.0075 

14.4 

9.0 

8.2 

Fhdo{Q) 

0.002 

18.0 

18.2 

18.7 

-F-NuII' (Q) 

0.00045 

27.0 

27.0 

27.0 

Th2o(Q) 

0.001 

16.0 

7.4 

oq 

06 


Table 3: cr parameters and R-factors for the different data sets in 
Case 1. 


3. Results and discussion 

3.1. Case 1: H/D isotopic substitution 

Four total scattering structure factors of four mixtures 
of light and heavy water were considered, as seen in Fig. 
2 of Ref. 0. The four compositions are: pure D 2 O, pure 



(T 

R 

AGD(r) 

0.000125 

57.7 

AGH(r) 

0.0005 

56.1 

AFd(Q) 

0.00003 

9.6 

AFh(Q) 

0.00002 

8.6 


Table 4: cr parameters and R-factors for the different data sets in 
Case 2. 

H 2 O, and two mixtures of light and heavy water with 
X = 0.5 (denoted as ’HDO’) and x = 0.64 (denoted as 
’Null’) {x is the mole fraction of light water in the mix¬ 
ture). In the case of the ’Null’ composition the hydrogen 
atoms have zero effective coherent scattering length and 
therefore, neutron scattering from this sample is sensitive 
to 0-0 correlations only. The partial radial distribution 
functions (goo, gon, and 5hh) were selected from Fig 4. 
of Ref. 0. Simulations were carried out to lit the 7 data 
sets simultaneously. 

The partial g{r) curves were fitted only in the inter¬ 
molecular regime (rmin = 2.24 A for goo, ?'min = 1-29 A 
for goH, and rmin = 1-89 A for 5 hh)- The total structure 
factors were fitted only from Qmin ~ 1-1 A“^, since the low 
Q region seemed to be the most uncertain (the Q range 
was the same for the four TSSF-s). (The quality of the 
fits over the low Q regime was also checked: it was found 
that below Qmin ~ 1-1 A”^ the curves fit poorly, so data 
from that region were omitted.) 

Two sets of the intramolecular O-H and H-H ranges 
(FNC-s) were checked (see Table[2]); note that these choices 
influence the molecular geometry that may form. The first 
one is stricter, corresponding to a realistic water molecule 
(Case la), whereas the second one (yielding ’elastic’ molecules. 
Case lb and Case Ic) leads to a wider range of intramolec¬ 
ular distances, in accord with the distance ranges depicted 
in Fig. 4 of Ref. 0 . 

The best fits obtained with realistic (as well as with 
’elastic’) molecules are shown in Fig. [TJ The goodness 
of fit values (R-factors) corresponding to the best fits are 
found in Table |31 The PRDF-s are displayed in Fig. [H 
Note that the high-Q parts of the H 2 O and D 2 O TSSF-s 
could not be approached sufficiently well when the ’strict’ 
definition of water molecules was in use. 

Comparing our gir) functions with Fig. 4 of Ref. 0 
it can be seen that in 0 the O-H and H-H intramolecular 
maxima are wider than those allowed by our ’strict’ FNC- 
s. Allowing such flexibility means that the molecules must 
suffer elastic deformations. In any case, we repeated our 
RMC calculations with these more flexible intramolecular 
constraints, see Case lb. Other simulation details (cutoffs, 
cr parameters, etc...) were the same as previously. The 
TSSF-s obtained by using ’elastic’ molecules are shown in 
Fig. [TJ and the corresponding R-factors are listed in Table 
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Q [Al 


Figure 1: RMC fits to the total scattering functions of the four water 
samples in Case 1. (Thin solid(blue) lines: RMC results for realistic 
molecules (Case la), dashed lines (red): RMC results for ’elastic’ 
molecules without ^^^^^(r) fit (Case lb), thick solid lines (green): 
RMC results for ’elastic’ molecules with (Case Ic), open 

circles (black) diffraction data from Fig. 2 of (^) 


[31 The R-factors are smaller, i.e. the fits are better when 
’elastic’ molecules are applied. However, the intramolecu¬ 
lar parts of the PRDF-s functions have become unrealis¬ 
tic, see Fig. |31 The H-H partial radial distribution func¬ 
tion has now two intramolecular peaks, at r « 1.38 and 
r R:! 1.65, and a local minimum at the expected distance 
(r Ri 1.55). The first peak is at the smallest value of the 
H-H intramolecular distance range allowed by the FNC-s 
(see Table HI). The deformation of the ’elastic’ molecule 
is also clearly visible in the distribution of the H-O-H (in¬ 
tramolecular) bond angles, see Fig. [3] for the realistic and 
’elastic’ molecules. The H-O-H angle distribution of real¬ 
istic molecules has a sharp peak around its most probable 
value (cca. 100°, cos0 ss -0.2). However, for the ’elastic’ 
molecule the H-O-H angle has a broad distribution with 
only a hump at around 103°, while the most developed 
(although still small) maximum is at 90°; this pattern is 
rather unrealistic. 

As a further effort to improve the situation, a third 
simulation (Case Ic) was also performed, by fitting also 
the intramolecular part of ^hh from Fig 4. of Q (note 
that the exact origin of this part remains undisclosed in 
Ref. 0). That is, in this RMC calculation 8 data sets 
were modeled simultaneously. The ’elastic’ definition of 
water molecules was used, while other simulation details 
(cutoffs, cr parameters) remained identical. The TSSF-s, 
PRDF-s and the distribution of H-O-H bond angles are 
presented in Figs (TJ H] and [31 R-factors are provided in 
Table |3l Clearly, the numerical values of the individual 
R-factors have increased somewhat, although there seem 
to be no visible difference in comparison with Case lb in 
terms of the TSSF-s and the intermolecular parts of the 
PRDF-s. The intramolecular part of 5 hh(?') now looks 
more sensible, although the H-O-H bond angle distribution 




Figure 2: Partial radial distribution functions obtained by RMC 
modeling. Thin solid (blue) lines: RMC results for realistic molecules 
(Case la), dashed lines (red): RMC results for ’elastic’ molecules 
without (Case lb), thick solid lines (green): RMC results 

for ’elastic’ molecules with (Case Ic), thick (black) lines 

with stars: RMC results for Case 2, open circles: intermolecular ^(r) 
functions from Fig. 4. of ^ open triangles: intramolecular Phh(^) 
function from Fig. 4. of Q. The insets show the intramolecular 
regions of the curves (b) 5 oh(^) Q^nd (c) 5 hh(?")- 




Figure 3: Distribution of the cosines of (a) H-O-H, (b) O- • • O- • • O, 
(c) H-O- • • O, (d) O-H- • • O angles. (Full squares: realistic molecules 
in Case la, open squares: ’elastic’ molecules without 
(Case lb), full triangle ’elastic’ molecules with fit (Case 

Ic), stars: Case 2. The lines are only guides to the eye.) 
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is still too wide to be realistic. 

It is instructive to notice that there are some Q regions 
where the experimental TSSF-s could not be approached 
as well for the ^H-containing samples (see particularly the 
’Null’ sample) as it is possible for pure D 2 O. This, together 
with the uncertainties of the intramolecular structure, may 
be considered as the manifestation of the difficult treat¬ 
ment of the incoherent inelastic background. 

For the characterization of the local structure in liq¬ 
uid water, the distributions of intermolecular O - ■ • O- ■ • O, 
H-0- ■ • O, and 0-H- • • O angles have also been calculated 
and shown in Fig. [3] for the three RMC simulations in 
conjunction with Case 1. A more detailed discussion will 
be provided after results concerning Case 2 are also in¬ 
troduced. Here it will suffice to note that although the 
molecular structure for these calculations (Cases la-c) are 
rather different, this does not seem to influence structural 
features connected to hydrogen bonding. 

3.2. Case 2: H/D and isotopic substitution 

The second set of data considered here is based on 
oxygen isotope substitution measurements 0, @)- Four 
different isotopic mixtures were investigated in these pub¬ 
lications: two heavy water (D 2^®0 and D 2 "’^*’ 0 ) and two 
light water (H2^®0 and H2''®'*0) samples. The first differ¬ 
ence functions from Fig. 4 and Fig. 5 of Ref. Q were 
taken for our calculations, both in Q and r space. These 
functions were calculated in the following way: after cor¬ 
rections the corresponding scattering cross sections were 
subtracted from each other, in order to eliminate the con¬ 
tributions from H-H or D-D correlations, together with 
the bulk of the inelastic scattering. The first difference 
functions are: 


AFuiQ) = F^^iQ) - = 

= Co(&i8o ~ &Lto) [SooiQ) ~ 1] + 

-I- 2cocd&d(& 18 Q — 6natQ )[5oD(g)-i] (1) 

and 

AFh(Q) = F^®(Q)-F“*(Q) = 

= Co(f)l 80 “ &nato) [SooiQ) ~ 1] + 

+ 2coch&h(&i8o ~ brn ^ to ) — 1] (2) 

where Cq and ba is the atomic fraction and bound co¬ 
herent neutron scattering length of chemical species a, 
SaisiQ) is a partial structure factor and Q = (47r/A)sin0 
is the magnitude of the scattering vector. Taking into 
account the neutron scattering lengths and the exact com¬ 
positions, the weighting factors for the “ 1] term 

in equations ([T]) and ([1]) were 0.00262 and 0.00263, and for 
the [5 'od(Q) — 1] and [5 'oh(Q) — 1] terms were 0.0059 and 
-0.0033, respectively. 



0 5 10 15 20 0 12,345 

Q [A’'] '' [A] 


Figure 4: RMC fits (black line) to the first difference functions of 
@1 (red dashed line), (a) AFj:,{r), (b) AFH(r)), (c) AGd(^), (d) 
AGh(^")- The insets show the region of the first peak in AGd('^) (c), 
and the first trough in AGh(’") (d). 

The corresponding real space functions AG'd(?') and 
AGh(?’) (which were obtained by Fourier transforming 
AFD(r) and AF'H(r)), appearing originally in Fig. 6 of 
0, were selected, too. These functions can be expressed 
from equations © and ([5]) as 

AGD(r) = 0.0059[goD(r) - f] + 0.00262[goo(r) - 1] (3) 

and 

AGH(r) =-0.0033[5oH(r-)-f]+0.00263[ffoo(r-)-l]. (4) 

RMC simulations were carried out with the aim of fit¬ 
ting the four data sets (as mentioned above) simultane¬ 
ously. Different starting configurations (random, as well 
as final configurations from Case f) were used. To find the 
best possible agreement between RMC and experiment, 
simulations were carried out in three different ways: (a) 
At first, only the two AF{Q) data sets were fitted, and 
the two AG(r)-s were taken into account in the second 
step, (b) The reverse way: AG(r)-s were applied in first 
step, (c) The four data sets were fitted together from the 
starting configurations. The best simultaneous fits for the 
four functions are shown in Fig. 01 while the correspond¬ 
ing partial radial distribution functions are given in Fig. 
m Goodness of fit (R-factor) values are in Table 01 The 
quality of the fits of the AG(r) data sets clearly shows 
that these four data sets are not entirely consistent: the 
largest discrepancies show up, interestingly, in terms of the 
’deuterated’ AGd(?')- 

Considering now the molecular structure, the situation 
is quite encouraging: the four target functions contain di¬ 
rect information on the intramolecular 0-H correlations 
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and clearly, the ’strict’ FNC values were easily applicable. 
The distribution of intramolecular H-O-H angles (see Fig. 
ED reflect far the most realistic molecular geometry: the 
maximum is symmetric, it is at the correct angle (above 
100 °), and there is no sign of any preference of either lim¬ 
iting value of the FNC range. (We note here that a variety 
of FNC-s have been tested here, too, with the hope that 
wider ranges would decrease differences between experi¬ 
mental and RMC results; as this desired outcome has not 
appeared, we have left the original ’strict’ boundaries un¬ 
changed.) 

Distributions of intermolecular O- • ■ O- • • O, H-0- ■ • O, 
and 0-H- • • O angles are also shown in Fig. El Each dis¬ 
tribution given for Case 2 is very similar to its Case 1 
counterpart: the largest discrepancies can be observed for 
the O- • • O- • • O angles where the maximum corresponding 
to angles of about 60° (close packing, hard sphere like fea¬ 
ture that is influenced by the nearest second neighbors) 
appears to be weakest for Case 2. This observation is con¬ 
sistent with differences between 0-0 PRDF-s (Fig. El) in 
the r range of 3.5 to 5 A. 

4. Conclusions 

Having completed a fairly detailed Reverse Monte Carlo 
modeling study of the two most recent, important sets of 
neutron diffraction data (^, 0 , the following findings are 
worth noting: 

(1) Concerning the data of Soper Q, both the TSSF-s 
and the intermolecular parts of the PRDF-s could be mod¬ 
eled satisfactorily, whereas the molecular structure had to 
be allowed to be somewhat unrealistically too v^ue. 

(2) Concerning the data of Zeidler et al. [1[ [3|i the 
level consistency between Q and r space information was 
somewhat inferior, but the molecular structure appeared 
as more realistic. 

(3) Clearly, even though the two sets of experiments 
(represented here by Case 1 and Case 2), as well as the 
two methods of data processing, could hardly be more 
different, essential details of the structure that concern 
hydrogen bonding turn out to be rather similar (in fact, 
nearly identical), as evidenced by Fig. El 

(4) The above finding provides some further hints as 
to why the structure of liquid water still presents difficul¬ 
ties: the duality (large differences between experimental 
data sets, and close similarities between conclusions on 
the microscopic structure; these appear at the same time) 
exposed by the present study leaves the interested reader 
somewhat puzzled. On the other hand, we may perhaps be 
assured that basic features of the hydrogen bonded struc¬ 
ture have been captured correctly over the years. 
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